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CASTRO: A New Compressible Astrophysical Solver. I. 
Hydrodynamics and Self- Gravity 
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M. J. Lijewski 1 , A. Nonaka 1 , M. Singer 2 , M. Zingale 4 

ABSTRACT 

We present a new code, CASTRO, that solves the multicomponent compress- 
ible hydrodynamic equations for astrophysical flows including self-gravity, nuclear 
reactions and radiation. CASTRO uses an Eulerian grid and incorporates adap- 
tive mesh refinement (AMR). Our approach to AMR uses a nested hierarchy of 
logically-rectangular grids with simultaneous refinement in both space and time. 



Q-c The radiation component of CASTRO will be described in detail in the next 

O ■ paper, Part II, of this series. 



Subject headings: methods: numerical, hydrodynamics, equation of state, gravi- 
tation, nuclear reactions 



Introduction 



q ■ In this paper, Part I of a two-part series, we present a new code, CASTRO, that solves 

the multicomponent compressible hydrodynamic equations with a general equation of state 
for astrophysical flows. Additional physics include self-gravity, nuclear reactions, and radi- 

•J~j ■ ation. CASTRO uses an Eulerian grid and incorporates adaptive mesh refinement (AMR). 

^ ■ Our approach to AMR uses a nested hierarchy of logically-rectangular grids with simulta- 

neous refinement of the grids in both space and time. Spherical (in ID), cylindrical (in ID 
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or 2D), and Cartesian (in ID, 2D or 3D) coordinate systems are supported. The radiation 
component of CASTRO will be described in detail in the next paper, Part II, of this series. 

There are a num ber of other adapti ve mesh c odes for compressib le astrophysi cal flows, 

most notably, ENZO (fo'Shea et a"Dl2005h . FLASH jFrvxell et al.lboooh . and RAGE JGittings et al 
20081 ). CASTRO differs from these codes in several ways. CASTRO uses an unsplit version of 
the piecewise parabolic method, PPM, with new limiters that avoid reducing the accuracy of 
the scheme at smooth extrema; the other codes are based on operator-split hydrodynamics, 
though the most recent release of FLASH, version 3.2, includes an unsplit MUSCL-Hancock 
scheme. The different methodologies also vary in their approach to adaptive mesh refine- 
ment. RAGE uses a cell- by-cell refinement strategy while the other codes use patch-based 
refinement. FLASH uses equal size patches whereas ENZO and CASTRO allow arbitrary 
sized patches. ENZO and FLASH enforce a strict parent-child relationship between patches; 
i.e., each refined patch is fully contained within a single parent patch; CASTRO requires 
only that the union of fine patches be contained within the union of coarser patches with a 
suitable proper nesting. Additionally, FLASH and RAGE use a single time step across all 
levels while CASTRO and ENZO support subcycling in time. All four codes include support 
for calculation of self-gravity 

It is worth noting that CASTRO us es the same grid stru cture as the low Mach number 
astrophysics code, MAESTRO (see, e.g., iNonaka et all (J2010h ) . This will enable us to map 
the results from a low Mach number simulation, such as that of the convective period and 
ignition of a Type la supernova, to the initial conditions for a compressible simulation such 
as that of the explosion itself, thus taking advantage of the accuracy and efficiency of each 
approach as appropriate. 



2. Hydrodynamics 

In CASTRO we evolve the fully compressible equations forward in time. The equations 
expressing conservation of mass, momentum, and total energy are: 



_ = -V • (pu) + S^p, 



d(pu) 

dt 
9{pE) 

dt 



-V • (puu) - Vj9 + pg + S cxt ,pu, 

-V • (puE + pu) + pH nnc + pu ■ g + S cxt , P E- 



(1) 
(2) 
(3) 



Here p, u, and E are the mass density, velocity vector, and total energy per unit mass, 
respectively. The total energy, E = e + u ■ u/2, where e is the specific internal energy. 
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The pressure, p, is defined by a user-supplied equation of state, and g is the gravitational 
acceleration vector. The source terms, S exttP , S ext>pu , and S cxt:P E are user-specified external 
source terms for the mass, momentum, and energy equations, respectively. For reacting 
flows, we evolve equations for mass fractions, X k : 

— 7j7 — = -V ■ (puX k ) + pu k + S cxttP x k . (4) 

where the production rates, u k , for species k are defined by a user-supplied reaction network. 
The reaction network also determines the energy generation rate H nuc . The mass fractions 
are subject to the constraint that ^2 k X k = 1. Again, a user-specified external source, 
S C xt, P x k , m ay be specified. Finally, CASTRO includes passively advected quantities, C£ dv , 
and auxiliary variables, C£ ux that satisfy 



d(pC? dv 



k I v7 ( /~iaAv\ 



V • (puCT) + S^cgt., (5) 

d(pCr) 



ot 



-V ■ (puCr) + S cxt , pC ^. (6) 



Advected and auxiliary variables are updated similarly, but they differ in their usage. In 
particular, auxiliary variables are passed into the equation of state routines. Examples 
of auxiliary and advected variables, respectively, might include the electron fraction, Y e , 
used in simulations of core collapse supernovae, and angular momentum in two-dimensional 
simulations of a rotating star in cylindrical (axisymmetric) coordinates. Both of these 
evolution equations include user-specified sources, S extpC ^d V and S cxttP c^ ux - We refer to 
U = (p, pu, pE, pX k , pC% dv , pC^) as the conserved variables. 



3. Equation of State and Reaction Network 

CASTRO is written in a modular fashion so that the routines for the equation of state 
and reaction network can be supplied by the user. However, for the test problems presented 
later we use routines that come with the CASTRO distribution. 

Each equation of state must provide an interface for obtaining thermodynamic quantities 
from p, e, and X k . One equation of state which is supplied with the CASTRO distribution 
is the gamma-law equation of state, which relates pressure and temperature, T, to p and e 
via: 

P = pe(l ~ 1) = — — ■ (7) 

pm p 

Here, 7, is the ratio of specific heats (e.g. 7 = 5/3 for a monatomic gas), ks is Boltzmann's 
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constant, m p is the mass of the proton, and the mean molecular weight, //, is determined by 



r?f< (8) 

with Ak the atomic weight of species k. 

The CASTRO distribution also includes more co mplex equations of state describing stel- 



lar m atter, including the Helmholtz equation of state (ITimmes fc Swestyll2000l ; lFryxeil et al. 



20001 ) which includes degenerate/relativistic electrons, ions (as a perfect gas), and radi- 
ation, and the Lattimer- Swesty equation of state, which describes dense nuclear matter 
(jLattimer fc Swestylll99ll ). For tabular equations of state, it is common that p, T, and X k 



are inputs, in which case a Newton-Raphson iteration is typically used to invert the equation 
of state. 

CASTRO can support any general reaction network that takes as inputs the density, 
temperature, and mass fractions, and returns updated mass fractions and the energy release 
(or decrease). The input temperature is computed from the equation of state before each 
call to the reaction network. In general, we expect the reaction network to evolve the species 
according to: 

^±=cu k (p,X k ,T). (9) 

at 

Reaction rates can be extremely temperature-sensitive, so in most cases, the reaction network 
should be written to evolve the temperature for the purposes of evaluating the rates. Close 
to nuclear statistical equilibrium, the energy release and change in abundances can rapidly 
change sign if the rate s are not evaluated with a temperature field consistent with the evolving 



energy (iMullerl 119861 ). At the end of the burning step, we use the energy release to update 



the total energy, E. The density remains unchanged during the burning. 



4. Gravity 

CASTRO supports several different options for how to specify and/or compute the 
gravitational acceleration. The simplest option is a gravitational field that is constant in 
space and time; this can be used for small-scale problems in which the variation of gravity 
throughout the computational domain is negligible. This option is available in ID Cartesian 
coordinates, 2D cylindrical or Cartesian coordinates, and 3D Cartesian coordinates. 

A second approach uses a monopole approximation to compute a radial gravity field 
consistent with the mass distribution. Because the algorithm subcycles in time we construct 
a separate ID radial density profile at each level at each time needed. Once the radial density 
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profile is defined, gravity is computed as a direct integral of the mass enclosed. This field is 
then interpolated back onto the original grids. 

The most general option is to solve the Poisson equation for self-gravity, i.e. solve 

V 2 = AirGp, (10) 

for </>, and define g = — V</>. This can be used in any of the coordinate systems. At bound- 
aries away from the star we set inhomogeneous Dirichlet boundary conditions for 0; these 
values are determined by computing the monopole approximation for g on the coarsest level, 
integrating this profile radially outward to create <j)(r), and interpolating onto the domain 
boundaries to define the boundary conditions for the solve. Boundaries that pass through 
the center of the star use symmetry boundary conditions. 

The Poisson equation is discretized using standard finite difference approximations and 
the resulting linear system is solved using geometric multigrid techniques, specifically V- 
cycles and red-black Gauss-Seidel relaxation. For multilevel calculations, special attention is 
paid to the synchronization of the gravitational forcing across levels, which will be discussed 
in Section [61 

There is also an option to add the gravitational forcing due to a specified point mass to 
either of the self-gravity options described above. 



5. Single-Level Integration Algorithm 

The time evolution of U can be written in the form 

_ = -V-F + S rC act + S, (11) 

where F is the flux vector, S rcact are the reaction source terms, and S are the non-reaction 
source terms, which includes any user-defined external sources, S cxt . We use Strang splitting 



(jStranglll968l ) to discretize the advection- reaction equations. In other words, to advance the 



solution, U, by one time step, At, we first advance the nuclear reaction network by At/2, 

U (1) = U" + ^SLct, (12a) 

then advect the solution by At, ignoring the reaction terms, 

U(2) = u(i) _ Atv . F n+y 2 + At 5 +° ? ( 12b ) 
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and finally advance the nuclear reaction network by another At/2, 



U 



n+l 



U(2) 



At. 



i(2) 



■"react • 



(12c) 



The construction of F is purely explicit, and based on an unsplit Godunov method. The 
solution, U, and source terms, S, are defined on cell centers; we predict the primitive vari- 
ables, Q = (p, u,p, pe, Xk, Cfc dv , Cfc UX ), from cell centers at time t n to edges at time t n+ b and 
use an approximate Riemann solver to construct fluxes, F™ +//2 , on cell faces. This algorithm 
is formally second-order in both space and time. 



5.1. Single-Level Flow Chart 



At the beginning of each time step, we assume that, in the case of self-gravity, g is 
defined consistently with the current mass distribution in U. The algorithm at a single level 
of refinement is composed of the following steps: 

Step 1: Advance the nuclear reaction network through a time interval of At/2. 

Define U^ = U n with the exception of 

At 



(pE) 



(i) 



(pXk) w 



(pE) n 

(px k y 



1" -^-(P-^nuc)", 



(13) 
(14) 



where (pH nuc ) n and (poJk) n are computed using calls to the user-defined reaction net- 
work. Note that p is unchanged during this step. 

Step 2: Advect the solution through At. 

Advance the solution using time-centered fluxes and an explicit representation of the 
source term, neglecting the contribution from reactions which are taken into account 
in Steps 1 and 4 (the asterisk superscript notation indicates that we will later correct 
this state to effectively time-center the source terms): 

U(2,**) = u(i) _ Atv . F n+y 2 + AtS (i)_ 



(15) 



where 



s(i) 
'u 



SpE 

s px k 

S pC *dv 

V Spc? u * ) 



(i) 



/ *->ext,p 

(pg) (1) + s cxt , pu 

(pu- g) (1) + S cxttP E 

Jcxt,pX k 

*-'cxt,pCJ dv 



(1) 



.s; 



(16) 



ext,pCg u 
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The construction of the fluxes is described in detail in Section 15.21 We note that in 
the single-level algorithm we can use the gravitational forcing computed in Step 3 of 
the previous time step, since the density has not changed. 

After the advective update, we ensure that the solution is physically meaningful by 
forcing the density to exceed a non-negative, user-defined minimum value. We also 
ensure that the mass fractions are all non-negative and sum to one. 

We also have an option for a user-defined sponge in order to prevent the velocities in 
the upper atmosphere from becoming too large, and subsequently, the time step from 
becoming too small. We multiply the velocity by 1/(1 + At k /dam P (p)), where k is the 
sponge strength, and /damp is a s mooth function of den sity that varies from to 1. Full 
details of the sponge are given in IZingale et al.l ( 120091 ). Finally, we adjust (pi?)^'**' to 
be consistent with u^ 2, **\ 

Step 3: Correct the solution with time- centered source terms and compute gravity at t n+l . 

We correct the solution by effectively time-centering the source terms. First, we correct 
U with updated external sources: 

At 



tt(2),* _ tt(2,**) i Ai ( e(2.**) c(i) \ 



(17) 



Next, we evaluate gravity using p( 2 '*\ If using full gravity we solve 

g(V) = _V0( 2 >*), V 2 ^ 2 '*) = 4nGp^*\ (18) 

where we supply an initial guess for </>( 2 '*) from the previous solve. In the single-level 
algorithm described here, g( 2 >*) is saved to be used as g^ in Step 2 of the next time 
step. This suffices in the single-level algorithm because p does not change between the 
end of Step 3 of one time step and the start of Step 2 of the next time step. 

We then correct the solution with the updated gravity: 

At 



(p£) (2) 



(P u 
(pE) 



,(2,*) 



+ 



(2,*) 



+ 



2 
At 



[ (pg )(2,*) _ (pg) (D] 
(pu-g) (2 '* } - (pu 



»(1) 



(19) 
(20) 



For all other conserved variables other than pu and pE^XJ^ = U*- 2 '*- 1 . We note here 
that the tim e discretization of the gravitatio nal forcing terms d iffers from that in 
the FLASH jFrvxell et all boooh and ENZO (fo'Shea et all l2005h codes, where the 
gravitati onal forcing at t n+ J 2 is computed by extrapolation from values at t n and t n_1 
(see also lBryan et al.lll995l ). Our discretization of the gravitational terms is consistent 
with our predictor-corrector approach in the handling of other source terms. 



Step 4: Advance the nuclear reaction network through a time interval of At/2. 
Define U n+1 = U^ with the exception of 



71+1 



(pE) 



(pX k ) n+1 



At 



(PE) (2) + -=-(pfli 



2 

At 



(pX k )W + —(pco 



»(2) 



I (2) 



(21) 
(22) 



We also include an option to modify any component of the new-time state as needed 
to account for special user requirements. 

Step 5: Compute the new time step. 

The time step is computed using the standard CFL condition for explicit methods, 
with additional constraints (such as one based on rate of burning) possible as needed. 
The user sets a CFL factor, cr CFL , between and 1. The sound speed, c, is computed 
by the equation of state, and for a calculation in n^im dimensions, 



At = a 



CFL 



min {At;} 

j=l...n dim 



where 



At; 



mm 



Ax 4 



|Uj| +c 
min x is the minimum taken over all computational grid cells in the domain. 



(23) 



(24) 



This concludes the single-level algorithm description. We note that whenever the kinetic 
energy dominates the total energy, making the calculation of e from E numerically unrelia ble, 
we use a method similar to the dual-energy approach described in iBryan et al.l (119951 ) to 
compute the internal energy with sufficient precision. In practice, this involves evolving pe 
in time and using this solution when appropriate. 



5.2. Construction of Fluxes 

W e use an unsplit Godun ov method with characteristic tracing and full corner coupling 
in 3D ( Miller fcColellal 120021) to compu te time-centered edge states. We have replaced the 
PPM limiters in lMiller &: Colellal (120021 ) with an updated PPM algorithm that is designed to 
preser ye accuracy at smooth extrema and is insensitive to a symmetries caused by roundoff 
error JColella fc Sekorall2008[ iMcCorquodale fc Colellalboioh. C ASTR O also has option s to 
use the unsplit piecewise-linear algorithm described in IColellal (Il990f ) ; ISaltzmanl (119941 ) , or 
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to retain the PPM limiters in iMiller Sz Colellal (120021 ) , which were originally developed in 
Colella fc Woodward! ( 119841 ) using a split integrator. 

There are four major steps in the construction of the face-centered fluxes, F n+ k that 
are used in Step 2 in Section 15.11 to update the solution. We also include details on the 
solution of the Riemann problem. In summary, 

Step 2.1: Rewrite the state, U^, in terms of primitive variables, Q*- 1 -*. 

Step 2.2: Construct a piecewise parabolic approximation of Q^ 1 * 1 within each cell. 

Step 2.3: Predict average values of Q^ on edges over the time step using charac- 
teristic extrapolation. 

Step 2.4: Compute fluxes, F n+ / 2 , using an approximate Riemann problem solver. 

We expand each of these steps in more detail below. 



Step 2.1: Compute primitive variables and source terms. 

We define Q^ = (p, u,p, pe,X k , C£ dv , C^ ux )^^. The pressure, p, is computed through 
a call to the equation of state using p, e, and X k . Note that we also include pe in 
Q; this quantity is used in the approximate Riemann solver to avoid an EOS call to 
evaluate the en e rgy fl ux, analogous to the effective dynamics for 7 = p/(pe) + 1 in the 



Colella &: Glazl ( 119851 ) approximate Riemann solver. 



For the overall integration algorithm, we want to include the effect of source terms 
except for reactions in the characteristic tracing (Step 2.2). (Reactions are treated 
separately using a symmetric operator split approach in Steps 1 and 4 of the algo- 
rithm.) The time evolution equations written in terms of the primitive variables, Q, 
and omitting contributions from reactions, are 



dQ 

dt 



( -u • Vp - pV • u \ 

-u • Vu - ±Vp 

— u • Vp — pc 2 ^J ■ u 

— u • V(pe) — (pe + p) V • u 

-u • VX fc 

-u ■ VCf v 



u ■ VC, 



aux 
k 



>Q 
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where 



s u 



>Q 



(-><: 



s, 
s 



O/^adv 






/ 



5*. 



cxt,p 

is 



~^cxt, pE + Pp^cxtp "I" ~~ ~ »Jcxt,pX fc 



cxt,pu 

P 



s, 
hs. 



cxt,pE 



V 



cxt,pX k 
p'-'cxt,pCJ dv 
-S e xt,pC% ux 



(25) 



7 



Here, c is the sound speed, defined as c = y/Tip/p, with Ti = dlogp/dlogp| s , with 
s the entropy. The remaining thermodynamic derivatives are p e = dp/de\ Pt x k , Vp = 
dp/dp\ et x k , and px k = dp/dXk\p, e ,Xj u ^ k) - Often, the equation of state is a function of 
p, T, and X&, and returns derivatives with these quantities held constant. In terms of 
the latter derivatives, our required thermodynamic derivatives are: 



Pe 



Pp 



Px k 




dp 



dX k 



P' T ' X j,U^k) 



de 



dp 



P,X k 



P,X k 

-1 



de 

dp 



T,X k 



dp 
df 



de 



cA 



dX k 



P' T ' X j,U^k) 



Step 2.2: Reconstruct parabolic profiles within each cell. 

In this step we construct a limited piecewise parabolic profile of each q in Q (we use 
q to denote an arbitrary primitive variable from from Q). These constructions are 
performed in each coordinate direction separately. The default option in CASTRO is 
to use a new limiting procedure that avoids reducing the order of the reconstruction at 
smoot h local extrema. The deta i ls of t his construction are given in IColella fc Sekora 
J2008I ): iMcCorquodale fc Colellal (fond ). In summary: 



Step 2.2a: For each cell, we compute the spatial interpolation of q n to the 
high and low faces of cell qi using a limited cubic interpolation formula. These 
interpolants are denoted by q i>+ and <&_. 

Step 2.2b: Construct quadratic profiles using qi-,qi, and q i+ . 



quad / \ 

Qi ( x ) 



%- + £{x) {ft,+ ~ ft,- + QaA 1 ~ £0*0]} 



(27) 
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<?e = Qqi - 3 



£(*) 



.r 



Hi 



h 



Qi,+ 



< £(x) < 1 



(28) 
(29) 



where h is the mesh spacing in the direction of interpolation. Also, as in Miller Sz Colella 
( 120021 ). we compute a flattening coefficient, x e [0)1]; used in the edge state 



prediction to further limit slopes near st rong shocks. The c omputation of x is 
identical to the approach used in FLASH f lFryxell et al.ll2000[ ). except that a flat- 
tening coefficient of 1 indicates that no additional limiting takes place, whereas a 
flattening coefficient of means we effectively drop order to a first-order Godunov 
scheme, which is opposite of the convention used in FLASH. 

Step 2.3: Characteristic extrapolation. 

We begin by extrapolating Q^ to edges at t n+ / 2 . The edge states are dual- valued, 
i.e., at each face, there is a left state and a right state estimate, denoted qL,i+V 2 an d 
(lR,i+ 1 / 2 (we write the equations in ID for simplicity). The spatial extrapolation is 
one-dimensional, i.e., transverse derivatives are omitted and accounted for later. 

• Step 2.3a: Integrate the quadratic profiles. We are essentially computing the 
average value swept out by the quadratic profile across the face assuming the pro- 
file is moving at a speed A&, where A^ is a standard wave speed associated with 
gas dynamics. 

Define the following integrals, where a^ = \Xk\At/h: 



2^,+ (o"fc 

Zi,-(0"fc 

Substituting (|2"T|) gives: 
Zi,+ (o~k) 



Y r(i-ya)h+<r k h 



qf ua (x)dx 
qf 1 ^ (x)dx 



(30a) 
(30b) 



0~k 



Qi,+ ~ Qi,- 

Qi,+ ~ Qi,- 



1 - g^fc ) qe,i 



1 - gO"fe ) q&,i 



(31a) 
(31b) 



Step 2.3b: Obtain a left and right edge state at £ n+ / 2 by applying a characteristic 
tracing operator (with flattening) to the integrated quadratic profiles. Note that 
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we also include the explicit source term contribution. 

[qi-I it+ (a k )]r k 



Ql^+Yz — Qi~Xi 2-^ h 

k:\ k >0 
k:\ k <0 



[qi -2j_(cr fc )]r* 



2 **' 

2 9,i ' 



(32a) 
(32b) 



In non-Cartesian coordinates, volume source terms are added to the traced states. 
Here, r^ and 1& are the standard right col umn and le ft row eigenvectors associated 
with the equations of gas dynamics (see iTorol 119971 ) . 

An unsplit approximation that includes full corner coupling is constructed by 
constructing increasingly accurate approximations to t he transverse derivatives . 
The details follow exactly as given in Section 4.2.1 in iMiller fc Colellal ( 120021 ). 
except for the solution of the Riemann problem, which is described in Step 2.4. 

Step 2.4: Compute fluxes 

The fluxes are computed using an appr oximate Riemann so lver. The solver used here 
is essentiall y the same as th at used in IColella et all (119971 ) , which is based on ideas 
discussed in iBell et all (119891 ) . This solver is computationally faster and considerabl y 
simpler than the approximate Riemann solver introduced by IColella fc Glaa (119851 ). 
The Colella and Glaz solver was based on an effective dynamics for 7 and was designed 
for real gases that are well-approximated by this type of model. The approximate 
Riemann solver used in CASTRO is suitable for a more general convex equation of 
state. 

As with other approximate Riemann solvers, an important design principle is to avoid 
additional evaluations of the equation of state when constructing the numerical flux. 
For that reason, we include pe in Q and compute (pe)L,R- The information carried in 
pe is overspecified but it allows us to compute an energy flux without an inverse call 
to the equation of state. 

The numerical flux computation is based on approximating the solution to the Riemann 
problem and evaluating the flux along the x/t = ray. The procedure is basically a 
two-step process in which we first approximate the solution in phase space and then 
interpret the phase space solution in real space. 

• Step 2.4a: To compute the phase space solution we first solve for p* and u*, the 
pressure between the two acoustic waves and the velocity of the contact discon- 
tinuity, respectively. These quantities are computed using a linearized approx- 
imation to the Rankine-Hugoniot relations. We first define I\ l/.r by using the 



From u* 


and p* 


we can 

* 
Pl,r 


compute 

p* 
= Pl,r + 


~PL,R 

2 '' 
C L,R 










(c1,r) 2 


= ^1,L,RP*l,r/ P*L,R-: 










(pe)*L,R 


= (pe)L,R + 


(P* ~ PL, 


R) 


+ p/p)l,r 

L L,R. 
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cell-centered values on either side of the interface. Next, we compute Lagrangian 
sound speeds, Wl = \/^\,lPlPl and W R = \/Ti >R p R p R and the corresponding 
Eulerian sound speeds c L>R = ^T ltLtR p LtR /p L>R . Then, 

* W L p R + W R p L + W L W R (u L - u R ) 
p = , (66a.) 

W L u L + W R u R + (p L - Pr) , , s 

" = wTTW R ' (33b) 



(34a) 

(34b) 

(34c) 

Vl,r v l,r, (34d) 

where v generically represents advected quantities (which includes transverse 
velocity components). Here, the notation L R refers to values on the left and right 
side of the contact discontinuity. 

Step 2.4b: The next step in the approximate Riemann solver is to interpret this 
phase space solution. If u* > then the contact discontinuity is moving to the 
right and numerical flux depends on the speed and structure of the acoustic wave 
connecting Q^ and Q* L associated with the A = u — c eigenvalue. Similarly if 
u* < then the contact is moving to the left and the numerical flux depends on 
the speed and structure of the acoustic wave connecting Q R and Q^ associated 
with the A = u + c eigenvalue. Here we discuss in detail the case in which u* > 0; 
the other case is treated analogously. 

For u* > we define \l = ul — cl and A^ = u* L — c* L . If p* L > pl then the wave is 
a shock wave and we define a shock speed a = ^(Xl + A£). For that case if a > 
then the shock is moving to the right and we define the Godunov state Qg = Ql', 
otherwise Qg = Q^- The rarefaction case is somewhat more complex. If both Al 
and X* L are negative, then the rarefaction fan is moving to the left and Qg = Q2- 
Similarly, if both A^ and A^ are positive, then the rarefaction fan is moving to the 
right and Qg = Ql- However, in the case in which \l < < X L , the rarefaction 
spans the x/t = ray and we need to interpolate the solution. For this case, we 
define 

Qg = «Q2 + (1-«)Q l (35) 
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where a = Al/(Al — X* L ). This choice of a corresponds to linearly interpolating 
Q through the rarefaction to approximate the state that propagates with zero 
speed. 

As noted above the case in which u* < is treated analogously. When u* = we 
compute Qg by averaging Q^ and Q R . For the Riemann problem approximation, 
we allow for user-specified floors for p, p and c to prevent the creation of non- 
physical values. 

The fluxes can then be evaluated from the final Qg- A small quadratic artificial 
viscosity that is proportional to the divergence of the velocity field is added to the flux 
in order to add additional dissipation at strong compressions. We also scale all the 
sp ecies fluxes so that they sum to the density flux, as in the sCMA algorithm described 
bv lPlewa fc Mtillerl fll999h . 



AMR 



Our approach to adaptive mesh refinement in CASTRO uses a nested hierarchy of 
logically-rectangular grids with simultaneous refinement of the grids in both space and time. 
The integration algorithm on the grid hierarchy is a recursive procedure in which coarse 
grids are advanced in time, fine grids are advanced multiple steps to reach the same time as 
the coarse grids and the data at different levels are then synchronized. 



The AMR methodology was introduced by iBerger fc Oligerl (119841); i t has been demon- 
strated to_beJughly successful for gas dynamics by lBerger fc Colellal (119891 ) in two dimensions 
and by iBell et all ( 119941 ) in three dimensions. 



6.1. Creating and Managing the Grid Hierarchy 

6.1.1. Overview 



The grid hierarchy is composed of different levels of refinement ranging from coarsest 
(£ = 0) to finest (£ = ^fi ncs t)- The maximum number of levels of refinement allowed, ^ max , is 
specified at the start of a calculation. At any given time in the calculation there may not 
be that many levels in the hierarchy, i.e. £fi ne st can change dynamically as the calculation 
proceeds as long as £fi ncs t < 4nax- Each level is represented by the union of non-overlapping 
rectangular grids of a given resolution. Each grid is composed of an even number of cells in 
each coordinate direction; cells are the same size in each coordinate direction but grids may 
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have different numbers of cells in each direction. Figure [T] shows a cartoon of AMR grids in 
two dimensions with two levels of refinement. 

In this implementation, the refinement ratio between levels £ and £ + 1, which we call 
r£, is always two or four, with the same factor of refinement in each coordinate direction. 
The grids are properly nested, in the sense that the union of grids at level £ + 1 is contained 
in the union of grids at level £. Furthermore, the containment is strict in the sense that, 
except at physical boundaries, the level £ grids are large enough to guarantee that there is a 
border at least n proper level £ cells wide surrounding each level £ + 1 grid (grids at all levels 
are allowed to extend to the physical boundaries so the proper nesting is not strict there). 
The parameter np roper is two for factor two refinement, and one for factor four refinement, 
since four ghost cells are needed for the PPM algorithm. 



6.1.2. Error Estimation and Regridding 



We initialize the grid hierarchy and regrid following the procedure outlined in lBell et al. 



( I1994I ). Given grids at level £ we use an error estimation procedure to tag cells where the 
error, as defined by user-specified routines, is above a given tolerance. Typical error criteria 
include first or second derivatives of the state variables or quantities derived from the state 
variables, or the state variables or derived quantities themselves. A user can specify that any 
or all of the criteria must be met to refine the cell; one can also specify criteria that ensure 
that a cell not be refined. For example, one could specify that a cell be refined if p > p cl -i t 
and ( (V 2 T) > (V 2 T) crit or |Vp| > |Vp| crit ), where p crit , (V 2 T) crit , and |Vp| crit are constants 
specified by the user. 

The tagged cells are grouped into rect angular grids at level £ using the clustering al- 



gorithm given in iBerger fc Rigoutsosl (Il99l[ ). These rectangular patches are refined to form 
the grids at level £ + 1. Large patches are broken into smaller patches for distribution to 
multiple processors based on a user-specified max_gridsize parameter. 

At the beginning of every ke level £ time steps, where kg > 1 is specified by the user at 
run-time, new grid patches are defined at all levels £ + 1 and higher if £ < £ max . In regions 
previously covered by fine grids the data is simply copied from old grids to new; in regions 
which are newly refined, data is interpolated from underlying coarser grids. 
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6.1.3. Enlarging the Domain 

The finest resolution of a calculation can vary in time; however, the coarsest resolution 
covering the domain does not change during a single run. However, a feature has been added 
to the CASTRO distribution that allows a user to restart a calculation in a larger domain 
covered by a coarser resolution, provided the data exists to initialize the larger domain. This 
is useful in simulations during which a star expands dramatically, for example. Using this 
strategy one could periodically stop the simulation, double the domain size, and restart the 
calculation in the larger domain. 



6.2. Multilevel Algorithm 

6.2.1. Overview 

The multilevel time stepping algorithm can most easily be thought of as a recursive 
procedure. In the case of zero or constant gravity, to advance level £, < £ < £ max the 
following steps are taken. Here the phrase, "Advance U" refers to Steps 1—4 of the single- 
level algorithm described in the previous section. 

• If £ = 0, compute the new time steps for all levels as follows 

— compute the appropriate time step for each level, At e ■* using the procedure de- 
scribed in Step 5 of the previous section, 

— define Rp as the ratio of the level cell size to the level £' cell size 

- define At = mxa v {R e At e '*), 

- define At 1 ' = Af/R v for all £' , < £' < £ max 

• Advance U at level £ in time as if it is the only level, filling boundary conditions for 
U from level £ — 1 if level £ > 0, and from the physical domain boundaries. 



If 



Advance U at level {£ + 1) for re time steps with time step At i+1 = —At 1 . 
Synchronize the data between levels £ and £ + 1 

* Volume average U at level £ + 1 onto level £ grids. 

* Correct U in all level £ cells adjacent to but not covered by the union 
of level £ + 1 grids thr ough an explicit refluxing operation as described in 



Berger fc Colellal (119891 ) 
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6.2.2. Monopole Gravity 

When we use the monopole gravity assumption in a multilevel simulation, we can no 
longer exploit the fact that p at level £ at the end of Step 3 of one time step is unchanged 
when one reaches the beginning of Step 2 of the next level £ time step. If £ < £ m axi then 
potential changes in p come from two sources: 

• p at level £ under the level £ + 1 grids is replaced by the volume average of p at level 

£ + 1; 

• the explicit refluxing step between levels £ and £ + 1 modifies p on all level £ cells 
adjacent to but not covered by the union of level £ + 1 grids. 

In addition, because the grids are dynamically created and destroyed through regridding, 
at the beginning of Step 2 of a level £ time step, there may not be a value for g from the 
previous step, because this region of space was previously not covered by level £ grids. 

In order to address all of these changes, we simply compute g^ at the beginning of 
Step 2 of each time step at each level, rather than copying it from g( 2 >*) from Step 3 of 
the previous time step as in the single-level algorithm. This captures any changes in grid 
structure due to regridding, and reflects any changes in density due to refluxing or volume 
averaging. 



6.2.3. Full Gravity Solve 

Overview 

Solving the Poisson equation for self-gravity on a multilevel grid hierarchy introduces ad- 
ditional complications. We start by defining some necessary notation. We define L e as an 
approximation to V 2 at level £, with the assumption that Dirichlet boundary conditions are 
supplied on the boundary of the union of level £ grids (we allow more general boundary 
conditions at physical boundaries), and define a level solve as the process of solving 

LV = 4tt<V 

at level £. 

We define Lf™ v as the composite grid approximation to V 2 on levels £ through m, and 
define a composite solve as the process of solving 

L c £ Z P 4> c ° mp = 47rGp comp 
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on levels £ through m. The solution to the composite solve satisfies 

L m comp = 4ivGp m 

at level m, but satisfies 

LV = 4ttG/ 

for £ < £' < m only on the regions of each level not covered by finer grids or adjacent to 
the boundary of the finer grid region. In regions of a level £' grid covered by level £' + 1 
grids the solution is defined as the volume average of the solution at H! + 1; in level £' cells 
immediately adjacent to the boundary of the union of level £' + 1 grid s, a modified in t erface 



operator is used that reflects the geometry of the interface (see, e.g., lAlmgren et al.l (Il998f ) 
for details of the multilevel cell-centered interface stencil). 

In an algorithm without sub cy cling; one can perform a composite solve at every time 



step, as described in iRickerl (120081 ) . to solve for on all levels. Because the CASTRO 
algorithm uses subcycling in time, however, we must use level solves at times when the 
solution is not defined at all levels, and then synchronize the solutions at different levels as 
appropriate. Even without changes in p due to volume averaging and refluxing, replacing a 
composite solve by separate level solves generates a mismatch in the normal gradient of at 
the boundary between each level. We correct these mismatches with a multilevel correction 
solve, which is a two-level composite solve for a correction to 0. In addition to correcting 
the solutions once the mismatch is detected, we add a correction term to later level solve 
solutions in order to minimize the magnitude of the correction that will be needed. 

Multilevel Algorithm 

At the start of a calculation, we perform a composite solve from level through £fi ncst to 
compute at all levels. In addition, after every regridding step that creates new grids at 
level £ + 1 and higher, a composite solve from level £ through £a nes t is used to compute at 
those levels. 



Following an approach similar to that described in iMiniati fc Colellal ( 120071 ). at the 



start and end of each level £ time step we perform a level solve to compute . The difference 
between 0™ mp and <j) e at the start of the time step is stored in 0^' corr . This difference is added 
to 0^ at the beginning and end of this level £ time step. Thus <p l + 0^' corr is identical to 0™ mp 
at the start of the time step; at the end of the time step it is an approximation to what the 
solution to the composite solve would be. In the event that the density does not change over 
the course of the time step, the effect of this lagged correction is to make £ + 0^ corr at the 
end of the time step identical to 0™ mp at that time, thus there is no mismatch between levels 
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to correct. In general, when the density is not constant, the effect of the lagged correction 
is to make the correction solve that follows the end of the time step much quicker. We now 
describe the two-level correction step. In the discussion below, we will refer to the two levels 
involved in a correction solve as the "coarse" and "fine" levels. 

At the end of r^ level t + 1 time steps, when the level £ + 1 solution has reached the 
same point in time as the level I solution, and after the volume averaging and refluxing steps 
above have been performed, we define two quantities on the coarse grid. The first is the 
cell-centered quantity, (<5p) c , which carries the change in density at the coarse level due only 
to refluxing. The second is the face-centered flux register, 

«i~* £+!>& <*> 

which accounts for the mismatch in the normal gradient of at coarse-fine interfaces. Here 
A c and A* represent area weighting factors on the coarse and fine levels, respectively. We 
define the composite residual, R comp , to be zero in all fine cells and in all coarse cells away 
from the union of fine grids, and 

R co ^ = ATrG(5p) c -(W-5F^\ c , (37) 

on all cells adjacent to the union of fine grids, where (V-)| c refers to the discrete divergence 
at the coarse level, where the only non-zero contribution comes from SF^ on the coarse-fine 
interface. We then solve 

L™-p 5<J) = R comp (38) 

and define the update to gravity at both levels, 

5g = -V(<50). (39) 

This update is used to correct the gravitational source terms. We define the new-time state 
after volume averaging but before refluxing as (p, u, pE, ...), and the contributions to the 
solution on the coarse grid from refluxing as ((5p) c , 5(pu) c , 5(pE) c , ...). Then we can define 
the sync sources for momentum on the coarse and fine levels, S^ c,c , and S s ^ c ^ , respectively 
as follows: 

^r ,c = (7f + (M c )(g c ' n+1 + ^g c )-p c g c ' n+1 
= [(M c g c,n+1 + (p c + (M c ) Sg c )} 

Sro/ = pf S gf. 
These momentum sources lead to the following energy sources: 

5 sync,c = S^-ivf + l/zAtcS^/r) 

s syn C J = S ^-{u-f + l/ 2 At f S^/f) 
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The state at the coarse and fine levels is then updated using: 

(pu) c ' n+1 = (pu) c + 5(pu) c + i/ 2 At c ^r' c , (pup +1 = (puY + i/aAt/S^'/, 
( P E) c ' n+1 = ( P E) C + 5{pE) c + y 2 At c ^ nc ' c , ( P E) f > n+1 = (p~E) f + y 2 At f S s / E ncJ . 

(The factor of 1/2 follows from the time-centering of the sources.) 
To complete the correction step, 

• we add 5(f) directly to <p e and <p e+1 and interpolate 5(f) to any finer levels and add it 
to the current at those levels. We note that at this point at levels £ and £ + 1 is 
identical to the solution that would have been computed using a two-level composite 
solve with the current values of density. Thus the new, corrected, at each level plays 
the role of comp in the next time step. 

• if level £ > 0, we transmit the effect of this change in to the coarser levels by updating 
the flux register between level £ and level £ — 1. In particular, we set 

SFS-^SF^ + ZA'Wgl. (40) 

Performance Issues 

The multilevel algorithm is not as computationally expensive as it might appear. Because 
multigrid is an iterative solver, the cost of each solve is proportional to the number of V- 
cycles, which is a function of the desired reduction in residual. We can reduce the number 
of V-cycles needed in two ways. First, we can supply a good initial guess for the solution; 
second, we can lower the desired reduction in residual. 

In the case of level solves, we always use from a previous level solve, when available, 
as a guess in the current level solve. Thus, even in a single-level calculation, from the 
beginning of the time step is used as a guess in the level solve at the end of the time step. 
If no regridding occurs, then at the end of one time step can be used as a guess for the 
level solve at the start of the next time step. The extent to which p changes in a time step 
dictates the extent to which a new computation of gravity is needed, but this also dictates 
the cost of the update. 

Similarly, there is no point in solving for 5(f) to greater accuracy than we solve for 0. 
When we do the correction solve for 5(f), we require only that the residual be reduced to 
the magnitude of the final residual from the level solve, not that we reduce the correction 
residual by the same factor. Thus, if the right-hand-side for the correction solve is already 
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small, the cost of the correction solve will be significantly less than that of the initial level 
solve. 



7. Software Design and Parallel Performance 

7.1. Overview 

CASTRO is implemented within the BoxLib framework, a hybrid C++ /Fortran90 
software system that provides support for the development of parallel structured-grid AMR 
applications. The basic parallelization strategy uses a hierarchical programming approach 
for multicore architectures based on both MPI and OpenMP. In the pure-MPI instantiation, 
at least one grid at each level is distributed to each core, and each core communicates with 
every other core using only MPI. In the hybrid approach, where on each socket there are 
n cores which all access the same memory, we can instead have one larger grid per socket, 
with the work associated with that grid distributed among the n cores using OpenMP. 

In BoxLib, memory management, flow control, parallel communications and I/O are 
expressed in the C++ portions of the program. The numerically intensive portions of the 
computation, including the multigrid solvers, are handled in Fortran90. The fundamental 
parallel abstraction in both the C++ and the Fortran90 is the MultiFab, which holds the 
data on the union of grids at a level. A MultiFab is composed of FAB's; each FAB is 
an array of data on a single grid. During each MultiFab operation the FAB's composing 
that MultiFab are distributed among the cores. MultiFab's at each level of refinement are 
distributed independently. The software supports two data distribution schemes, as well as a 
dynamic switching scheme that decides which approach to use based on the number of grids 
at a level and the numbe r of processors. Th e first scheme is based on a heu ristic knapsack 



algorithm as described in ICrutchfieldl ( 119911 ) and in lRendleman et al.1 (120001 ). The second is 



based on the use of a Morton-ordering space- filling curve. 

Each processor contains meta-data that is needed to fully specify the geometry and 
processor assignments of the MultiFab's. At a minimum, this requires the storage of an 
array of boxes specifying the index space region for each AMR level of refinement. One of 
the advantages of computing with fewer, larger grids in the hybrid OpenMP-MPI approach 
is that the size of the meta-data is substantially reduced. 
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7.2. Parallel Output 



Data for checkpoints and analysis are written in a self-describing format that consists of 
a directory for each time step written. Checkpoint directories contain all necessary data to 
restart the calculation from that time step. Plotfile directories contain data for postprocess- 
ing, visualization, and analytics, which can be read using amrvis, a cust omized visualization 



packa ge developed at LBNL for visualizing data on AMR grids, or Visit (j Visit User's Manual 



20051 ). Within each checkpoint or plotfile directory is an ASCII header file and subdirectories 
for each AMR level. The header describes the AMR hierarchy, including number of levels, 
the grid boxes at each level, the problem size, refinement ratio between levels, step time, 
etc. Within each level directory are the MultiFab files for each AMR level. Checkpoint and 
plotfile directories are written at user-specified intervals. 

For output, each processor writes its own data to the appropriate MultiFab files. The 
output streams are coordinated to only allow one processor to write to a file at one time 
and to try to maintain maximum performance by keeping the number of open data streams, 
which is set at run time, equal to the number of files being written. Data files typically 
contain data from multiple processors, so each processor writes data from its associated 
grid(s) to one file, then another processor can write data from its associated grid(s) to that 
file. A designated I/O Processor writes the header files and coordinates which processors are 
allowed to write to which files and when. The only communication between processors is for 
signaling when processors can start writing and for the exchange of header information. We 
also use the C++ setbuf function for good single file performance. While I/O performance 
even during a single run can be erratic, recent timings on the Franklin machine (XT4) at 
NERSC indicate that CASTRO's I/O performance, when run with a single level composed 
of multiple uniform ly-sized grids, matches some of the top re sults for the N5 IOR benchmark 



roughly 13GB/s) (IFranklin Performance Monitoring! |2010| ). For more realistic simulations 



with multiple grids at multiple levels, CASTRO is able to write data at approximately 5 
GB/s sustained, over half of the average I/O benchmark reported speed. 



7.3. Parallel Restart 

Restarting a calculation can present some difficult issues for reading data efficiently In 
the worst case, all processors would need data from all files. If multiple processors try to read 
from the same file at the same time, performance problems can result, with extreme cases 
causing file system thrashing. Since the number of files is generally not equal to the number 
of processors and each processor may need data from multiple files, input during restart is 
coordinated to efficiently read the data. Each data file is only opened by one processor at 
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a time. The IOProcessor creates a database for mapping files to processors, coordinates the 
read queues, and interleaves reading its own data. Each processor reads all data it needs 
from the file it currently has open. The code tries to maintain the number of input streams 
to be equal to the number of files at all times. 

Checkpoint and plotfiles are portable to machines with a different byte ordering and 
precision from the machine that wrote the files. Byte order and precision translations are 
done automatically, if required, when the data is read. 



7.4. Parallel Performance 

In Figure [2] we show the scaling behavior of the CASTRO code, using only MPI- 
based parallelism, on the jaguarpf machine at the Oak Ridge Leadership Computing Facility 
(OLCF). A weak scaling study was performed, so that for each run there was exactly one 
64 3 grid per processor. We ran the code with gravity turned off, with the monopole approx- 
imation to gravity, and with the Poisson solve for gravity. The monopole approximation to 
gravity adds very little to the run time of the code; with and without the monopole approxi- 
mation the code scales excellently from 8 to 64,000 processors. For the 64,000 processor case 
without gravity, the time for a single core to advance one cell for one time step is 24.8 fis. 

Good scaling of linear solves is known to be much more difficult to achieve; we report 
relatively good scaling up to only 13,824 processors in the pure-MPI approach. An early 
strong scaling study contrasting the pure-MPI and the hybrid-MPI-OpenMP approaches 
for a 768 3 domain shows that one can achieve at least a factor of 3 improvement in linear 
solver time by using the hybrid approach at large numbers of processors. Improving the 
performance of the linear solves on the new multicore architectures is an area of active 
research; more extensive development and testing is underway. 

We also ran a scaling study with a single level of local refinement using the monopole 
gravity approximation. In this MPI-only study, there is one 64 3 grid at each level for each 
processor. Because of subcycling in time, a coarse time step consists of a single step on 
the coarse grid and two steps on the fine grid. Thus, we would expect that the time to 
advance the multilevel solution by one coarse time step would be a factor of three greater 
than the time to advance the single-level coarse solution by one coarse time step, plus any 
additional overhead associated with AMR. From the data in the figure we conclude that 
AMR introduces a modest overhead, ranging from approximately 5% for the 8 processor 
case to 19% for the 64,000 processor case. By contrast, advancing a single-level calculation 
at the finer resolution by the same total time, i.e., two fine time steps, would require a factor 
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of 16 more resources than advancing the coarse single-level solution. 



8. Test Problems 

In this section we present a series of calculations demonstrating the behavior of the hy- 
drodynamics, self-gravity, and reaction components of CASTRO. The first set contains three 
one-dimensional shock tube problems, including Sod's problem, a double rarefaction prob- 
lem, and a strong shock problem. We follow this with Sedov- Taylor blast waves computed in 
ID spherical coordinates, 2D cylindrical and Cartesian coordinates, and 3D Cartesian coor- 
dinates. Our final pure-hydrodynamics test is a 2D Rayleigh- Taylor instability. We use this 
problem to contrast the differences in the flow found using dimensionally split and unsplit 
methods with piecewise linear, PPM with the old limiters, and PPM with the new limiters. 

We then present two examples that test the interaction of the self-gravity solvers with 
the hydrodynamics in 3D Cartesian coordinates. In the first case a star is initialized in 
hydrostatic equilibrium and we monitor the maximum velocities that develop; in the sec- 
ond, the homologous dust collapse test problem, a uniform- density sphere is initialized at 
a constant low pressure, and collapses under its own self-gravity These tests more closely 
examine the 3D spherical behavior we expect to be present in simulations of Type la and 
Type II supernovae. 

We perform a test of the coupling of the hydrodynamics to reactions. This test consists 
of a set of buoyant reacting bubbles in a stratified stellar atmosphere. We compare the 
CASTRO results to those of the FLASH code. 

Finally, we note that a previous comparison of CASTRO t o our low Mach number 



hydrodynamics code, MAESTRO, can be found in iNonaka et al.l (120101 ). In that test, we 
took a 1-d spherical, self-gravitating stellar model and watched it hydrostatically adjust as 
we dumped energy into the center of the star. The resulting temperature, pressure, and 
density profiles agreed very well between the two codes. 



8.1. Shock Tube Problems 

To test the behavior of the hydrodynamics solver, we run several different ID shock tube 
problems. The setup for these problems consists of a left and right state, with the interface 
in the center of the domain. All calculations use a gamma- law equation of state with 7 = 1.4. 
We show results from each problem run using ID Cartesian coordinates, but we have verified 
that the results are identical when each problem is run in 2D or 3D Cartesian coordinates 
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and the interface is normal to a coordinate axis. The length of the domain is always taken 
as 1.0, with the interface in the center. We use a base grid of 32 cells, with two additional 
levels of factor 2 refinement, for an effective resolution of 128 cells. The refinement criteria 
are based on gradients of density and velocity. In the case of the double rarefaction we 
also present results from runs with two levels of factor 4 refinement (effective resolution of 
512 cells) and three levels of factor 4 refinement (effective resolution of 2 048 ce l ls). In each 
case, analytic solutions are found using the exact Riemann solver from iTorol (119971 ). All 
calculations are run with the new PPM limiters and a CFL number of 0.9. For each problem 
we show density, pressure, velocity, and internal energy. 



8.1.1. Sod's Problem 



The Sod problem (jSodl Il978l ) is a simple shock tube problem that exhibits a shock, 
contact discontinuity, and a rarefaction wave. The non-dimensionalized initial conditions 
are: 

PL = 1 PR = 0.125 

u L = u R = (41) 

PL = 1 PR = 0.1 

This results in a rightward moving shock and contact discontinuity, and a leftward moving 
rarefaction wave. Figure |3] shows the resulting pressure, density, velocity, and internal energy 
at t — 0.2 s. We see excellent agreement with the exact solution. 



8.1.2. Double Rarefaction 



The double rarefaction problem tests the behavior of the hydro dyn amics algorit hm in 
regions where a vacuum is created. We run the problem as described in ITorol (119971 ). The 
non-dimensionalized initial conditions are: 



PL 

I'L 
PL 



1 

-2 
0.4 



PR = 1 

u R = 2 
Pr = 0.4 



(42) 



This results in two rarefaction waves propagating in opposite directions away from the center. 
As a result, matter is evacuated from the center, leaving behind a vacuum. Figure H] shows 
the CASTRO solutions at t — 0.15 s. The agreement with the exact solution is excellent at 
the 128-cell resolution for density, pressure and velocity; the internal energy is more sensitive, 
but clearly converges to the analytic solution except at the center line. This is a very common 
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pathology for this problem, since the internal energy, e, is derived from equati on ([71) usin g 



values of p and p which are both approaching zero in the center of the domain fjToro 



19971 ). 



8.1.3. Strong Shock 

The fin al shock-tub e problem we try is a strong shock. We initialize the problem as 
described in iTorol (J1997I ). The initial conditions are: 



PL = 1 

u L = 
p L = 1000 



PR = 1 

u R = 
PR = 0.01 



(43) 



The initial pressure jump of six orders of magnitude results in a strong rightward moving 
shock. This large dynamic range can cause trouble for some hydrodynamics solvers. The 
shock is followed very closely by a contact discontinuity. A leftward moving rarefaction is 
also present. Figure [5] shows the CASTRO results at t — 0.012 s. We see good agreement 
between the CASTRO results and the exact solution. 



8.2. Sedov 

Another standard hydrodynamics test is the Sedov- Taylor blast wave. The problem 
setup is very simple: a large amount of energy is deposited into the center of a uniform 
domain. This drives a blast wave (spherica l or cy li ndric al, depending on the domain geom- 
etry). An a nalytic solution is provide d by ISedovl ( 119591 ). We use a publicly available code 
described by lKamm fc Timmesl (120071 ) to generate the exact solutions. 



The Sedov explosion can test the geometrical factors in the hydrodynamics scheme. A 
cylindrical blast wave (e.g. a point explosion in a 2D plane) can be modeled in 2D Cartesian 
coordinates. A spherical blast wave can be modeled in ID spherical, 2D axisymmetric 
(cylindrical r-z), or 3D Cartesian coordinates. 

In the Sedov problem, the explosion energy, S exp (in units of energy, not energy/mass 
or energy/ volume), is deposited into a single point, in a medium of uniform ambient density, 
Pambient; an d pressure, Pambicnt- Initializing the problem can be difficult because the small 
volume is typically only one cell in e xtent, which can lead to grid imprinting in the solution. A 
standard approach (see for example iFryxell et al.ll2000l : lQmang et al.ll2006l and the references 
therein) is to convert the explosion energy into a pressure contained within a certain volume, 
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V init , of radius r init as 

PMl = <lZ^E. (44) 

•'init 

This pressure is then initialized to p init in all of the cells where r < r init . We use the 
gamma- law equation of state with 7 = 1.4. 

To further minimize any grid effects, we do subsampling in each cell: each cell is divided 
it into iV su b subcells in each coordinate direction, each subcell is initialized independently, and 
then the subcells are averaged together (using volume weighting for spherical or cylindrical 
coordinates) to determine the initial state of the full cell. 

For these runs, we use p am bient = 1 g cm -3 , p am bicnt = 10~ 5 dyn cm -2 , £ exp = 1 erg, 
rinit = 0.01 cm, and iV su b = 10. A base grid with Ax = 0.03125 cm is used with three levels 
of factor 2 refinement. For most geometries, we model the explosion in a domain ranging 
from to 1 cm in each coordinate direction. In this case, the base grid would have 32 
cells in each coordinate direction and the finest mesh would correspond to 256 cells in each 
coordinate direction. For the 2D axisymmetric case, we model only one quadrant, and the 
domain ranges from to 0.5 cm. All calculations were run with a CFL number of 0.5, and 
the initial time step was shrunk by a factor of 100 to allow the point explosion to develop. 
We refine on regions where p > 3 g cm -3 , Vp > 0.01 g cm -3 cm -1 , p > 3 dyn cm -2 , or 
Vp > 0.01 dyn cm -2 cm -1 . 

Figure |6] shows the CASTRO solution to a spherical Sedov explosion at time t = 0.01s, 
run in ID spherical, 2D cylindrical, and 3D Cartesian coordinates. For the 2D and 3D 
solutions, we compute the radial profile by mapping each cell into its corresponding radial 
bin and averaging. The radial bin width was picked to match the width of a cell at the 
finest level of refinement in the CASTRO solution. The density, velocity, and pressure plots 
match the exact solution well. As with the double rarefaction problem, the internal energy 
is again the most difficult quantity to match due to the vacuum region created at the origin. 
Figure [7] shows the same set of calculations run with 4 levels of factor 2 refinement. Here 
the agreement is even better. Figure |8] shows the CASTRO solution at time t — 0.1s to a 
cylindrical Sedov explosion, run in 2D Cartesian coordinates. 



8.3. Rayleigh- Taylor 

The Rayleigh- Taylor i nstability results when a dense fluid i s placed over a less-dense 



fluid in a gravitational field ( JTaylorlll950l ; lLayzerlll955l ; ISharplll984l ). The interface is unstable 



and a small perturbation will result in the growth a buoyant uprising bubbles and dense, 
falling spikes of fluid. This instability provides a mechanism for mixing in many astrophysical 
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systems. Despite its seemingly simplistic nature, only the linear growth regime is understood 
analytically (see for example IChandrasekharl Il96ll ). In the non-linear reg ime, Rayleigh- 
Taylo r instability calculations are often used as a means of code validation ( IDimonte et al. 
2004h . 



For our purposes, the R-T instability provides a good basis to compare different choices 
of the advection algorithm. We model a single-mode Rayleigh- Taylor instability — a pertur- 
bation consisting of a single wavelength that disturbs the initial interface. Short-wavelength 
perturbations have a faster growth rate than long- wavelength perturbations, so grid effects 
can easily drive the instability on smaller scales than our initial perturbation. No viscous 
terms are explicitly modeled. 

We choose the density of the dense fluid to be p 2 = 2 g cm" 3 and the light fluid is 
Pi = 1 g cm -3 . The gravitational acceleration is taken to be g = — 1 cm s -2 in the vertical 
direction. The gamma-law equation of state is used with 7 = 1.4. Our domain has a width 
of L x = 0.5 cm and a height of L y = 1 cm. The initial interface separating the high and low 
density fluid is centered vertically at L y /2, with the density in the top half taken to be p 2 
and the density in the lower half p\. Since g and pi, p 2 are constant, we can analytically 
integrate the equation of hydrostatic equilibrium to get the pressure in both the high and 
low-density regions of the domain: 



p(y) 



Pbasc + pigy y < L y /2 

Pbase + PigLy/2 + p 2 g(y - L y /2) y > L y /2 



(45) 



where y is the vertical coordinate, and pbase is the pressure at the base of the domain. We 
take p baS e = 5 dyn cm -2 . 

To initiate the instability, the interface is perturbed by slightly shifting the density, 
keeping the interface centered vertically in the domain. We define the perturbed interface 
height, tp, to be a function of position in the x-direction as 



ip(x) 



A 
~2 



2irx\ f2ir(L x — x] 

I -= — ) + cos 



(46) 



with the amplitude, A = 0.01 cm. We note that the cosine part of the perturbation is done 
symmetrically, to prevent roundoff error from introducing an asymmetry in the flow. The 
density is then perturbed as: 



P(x,y) = pi + 



P2-P1 



tanh 



h 



(47) 



The tanh profile provides a slight smearing of the initial interface, over a smoothing length 
h. We take h = 0.005 cm. 
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In Figure [H we present simulation results for the Rayleigh- Taylor problem at t — 2.5s 
for several different variants of the hydrodynamics. All calculations were run with 256 x 512 
grid cells. In the bottom right image we show the results obtained using the unsplit PPM 
with the new limiter used in CASTRO. The left and middle images on the bottom row 
are results using the u nsplit piecewise linear method and unsplit PPM with limiters as in 



Miller fc Colellal (120021 ). respectively. The results with all three methods are reasonably 
good; however, the piecewise linear and original PPM limiter both exhibit mild anomalies 
at the tip of both the bubble and the spike. 

In the upper row, we present results for the Rayleigh- Taylor problem using operator- 
split analogs of the unsplit methods. The details of the algorithms such as limiters, Riemann 
solver, etc. are the same as in the unsplit methods; the only difference is the use of operator 
splitting. We note that all three of the operator-split methods produce spurious secondary 
instabilities. This behavior is a direct result of the operator-split approach. Physically, 
for these low Mach number flows, the density field is advected by a nearly incompressible 
flow field, and remains essentially unchanged along Lagrangian trajectories. However, in 
regions where there is significant variation in the local strain rate, an operator-split inte- 
gration approach alternately compresses and expands the fluid between subsequent sweeps. 
This alternating compression / expansion provides the seed for the anomalies observed with 
operator-split methods. 

We note that both the CPU time and the memory usage are roughly a factor of two larger 
for the unsplit algorithm than for the split algorithm in this two-dimensional implementation. 
For a pure hydrodynamics problem with gamma-law equation of state this factor is nontrivial; 
for a simulation that uses the full self-gravity solver, a realistic reaction network, a costly 
equation of state, or significant additional physics, the additional cost of the hydrodynamic 
solver may be negligible. 

In 3D one might expect the ratio of CPU time for the unsplit algorithm relative to 
the split algorithm to be be even larger than in 2D because of the additional Riemann 
solves required to construct the transverse terms. However, this effect is counterbalanced 
by the need to advance ghost cells in the split algorithm to provide boundary conditions 
for subsequent sweeps. Consequently, we observe an increase in CPU time that is slightly 
less than the factor of two observed in 2D. The 3D implementation of the unsplit algorithm 
in CASTRO uses a strip-mining approach that only stores extra data on a few planes at a 
time, so we see an increase of less than 10% in the memory required for the unsplit integrator 
compared to the split integrator in 3D. 
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8.4. Stationary Star Gravity 

A challenging problem for a hydrodynamics code is to keep a star in hydrostatic equi- 
librium. Because of the different treatment of the pressure, density, and gravitational ac- 
celeration by the hydrodynamics algorithm, small motions can be driven by the inexact 
cancellation of Vp and pg. This is further exaggerated by modeling a spherical star on a 3D 
Cartesian grid. Here we test the ability of CASTRO to maintain hydrostatic equilibrium for 
a spherical, self-gravitating star. 

Our initial model is a nearly-Chandrasekhar mass, carbon-oxygen white dwarf, which 
is generated by specifying a core density (2.6 x 10 9 g cm -3 ), temperature (6 x 10 8 K), 
and a uniform composition (A( 12 C) = 0.3,AT( 16 O) = 0.7) and integrating the equation of 
hydrostatic equilibrium outward while constraining the specific entropy, s, to be constant. 
In discrete form, we solve: 

Poj+i = Po,j + 2 Ar (Po,j + Po,j+i)9j+y 2 , (48) 

s o,j+i — s o,j, (49) 

with Ar = 1.653125 x 10 5 cm. We begin with a guess of Poj+i an d ^oj+i an d use the equation 
of state and Newton-Raphson iterations to find the values that satisfy our system. Since 
this is a spherical, self-gravitating star, the gravitation acceleration, <7j+y 2 , is updated each 
iteration based on the current value of the density. Once the temperature falls below 10 7 K, 
we keep the temperature constant, and continue determining the density via hydrostatic 
equilibrium until the density falls to 10~ 4 g cm" 3 , after which we hold the density constant. 
This uniquely determines the initial model. We note that this is the same procedure we follow 
to initialize a connecting white dw arf for the multilevel low Mach number code, MAESTRO, 



described in iNonaka et al.l (J2010f l. 

We map the model onto a (5 x 10 8 cm) 3 domain with 192 3 , 384 3 , and 768 3 grid cells, 
and center the star in the domain. We let the simulation run to 1 s, and compare the 
maximum magnitude of velocity vs. time and the magnitude of velocity vs. radius at t = 1 s, 
a time greater than two sound-crossing times. We only consider regions of the star at 
r < 1.8 x 10 8 cm, which corresponds to a density of p w 5.4 x 10 5 g cm -3 . Note that the 
density reaches the floor of 10~ 4 g cm -3 at r = 1.9 x 10 8 cm. We turn on the sponge at the 
radius where p = 100 g cm -3 and the sponge reaches its full strength at the radius where 
p = 10~ 4 g cm -3 with a sponge strength of k = 10 00 s -1 . We use a CFL of 0.9 and no 



refine ment. We use the Helmholtz equation of state (jTimmes fc Swestyl |2000| ; iFryxell et al. 



20001 ) and no reactions are modeled. 

Figure [10] shows a plot of the maximum magnitude of velocity vs. time. At each of the 
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three resolutions, we show the results using a monopole gravity approximation and Poisson 
solve for gravity. We note that in each simulation, the maximum velocity is not strictly 
increasing, leading us to believe that over longer periods of time the velocities will remain 
small. We note that sound speed at the center of the star is approximately 9.4 x 10 8 cm/s, so 
at the highest resolution, the peak velocity is less than 1% of the sound speed. The monopole 
and Poisson cases match up very well, except for the finest resolution. The reason why we 
see larger peak velocities in the finest resolution Poisson solver simulation is due to the large 
velocities at the edge of the star. 

Figure [11] shows a plot of the magnitude of velocity vs. radius at t — 1 s. Again, at each 
of the three resolutions, we show the results using a monopole gravity approximation and 
Poisson solve for gravity. Here, we see clear second order convergence in the max norm, and 
the monopole and Poisson simulations agree best at the highest resolution. We also see how 
in the finest resolution runs, the velocities at the edge of the star can become large, but this 
is likely outside the region of interest for a typical simulation. 



8.5. Homologous Dust Collapse 

As a second test of the gravity solver in CASTRO we implement the homologous dust 
collapse test problem, a 'pressure-less' configuration that collapses under its own self-gravity 
A n analytic solution that describes the radius of the sphere as a function of time is found 



in 



Colgate fc White! (11966). Ou r implementation of this probl em follows that described in 



FLASH 3.2 User's Guide! (120091 ): Monchmever &: Mulled (119891 ). The problem is initialized 
with a sphere with a large, uniform density, po> of radius r . The pressure everywhere should 
be negligible, i.e., the s ound crossing time should be m u ch longer than th e free-f all collapse 
time (see, for example, IfLASH 3.2 User's Guidell2009l ). JColgate fc White! jl966h use p = 0. 
We choose a value that does not appear to affect the dynamics. As the sphere collapses, the 
density inside should remain spatially constant, but increase in value with time. 



10 9 g cm 3 and Tq 



Following IFLASH 3.2 User's Guide! (1200a ). we take p 
10 8 cm. The pressure is not specified, so we take it to be 10 15 dyn cm -2 



6.5 x 
Outside of the 



sphere, we set the density to p a mbicnt = 10 5 g cm 3 . Finally, since the sharp cutoff at the 
edge of the sphere is unphysical, we smooth the initial profile by setting 



P = Po 



P0 — Pambient 



1 + tanh 



r — tq 

h 



(50) 



with the smoothing length, h = 4 x 10 6 <C ro. We use the gamma-law equation of state with 
7 = 1.66. 
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Figure 1121 shows the radius vs. time for the ID, 2D, and 3D simulations as compared to 
the exact solution. In all three cases we see excellent agreement with the exact solution. 



8.6. Reacting Bubbles in a Stellar Atmosphere 

A final test is a code comparison of the evolution of three reacting bubbles in a plane- 
paralle l stellar atmosp h ere. T his problem is almost identical to the setup described in Section 
4.2 of lAlmgren et al.l (120081 ) with two minor differences. First, we eliminate the stably 
stratified layer at the base of the atmosphere by setting the lower y extrema of the domain 
to 5.00625 x 10 7 cm — this way, the bottommost row of cells in the domain is initialized 
with the specified base density (2.6 x 10 9 g cm -3 ) and temperature. Second, we set the base 
temperature of the atmosphere to 6 x 10 8 K (instead of 7 x 10 8 K) to minimize the amount of 
reactions occurring near the lower domain boundary. Three temperature perturbations are 
seeded in pressure- equilibrium wit h a ran ge of heights and widths as specified by equation 
(87) and Table 1 of lAlmgren et al.l (120081 ). We use a uniform computation grid of 384 x 576 
cells and a domain width of 2.16 x 10 8 cm. 



We compare the evolution to the FLASH code (JFryxell et al.ll2000l ). version 2.5, using 
the standard dimensionally-split PPM hydrodynamics module that comes with FLASH. 
The lower boundary condition in both cases provides hydrostatic support by integrating the 
equation of hydrostatic equilibrium together with t he equations of state into the ghost cells, 
assuming a constant temperature, as described in IZingale et al.l (120021 ). The left and right 
boundary is periodic. We u s e the same single step ( 12 C + 12 C — > 24 Mg) reaction module 
described in lAlmgren et al.l (120081). Both codes use the g eneral stellar equation of state 



described in iFryxell et al.l (120001 ): iTimmes fc Swestyl (120001 ) with the Coulomb corrections 
enabled. 



Figures [TBI and [T41 show contours of the temperature and AT( 24 Mg) after 2.5 s of evolution 
for both FLASH and CASTRO. We see excellent agreement between the two codes in terms 
of bubble heights and contour levels. 



8.7. Type la Supernova 

As a final example, in Figure [151 we show a 2D snapshot of tempe rature from a 3D 
calculation of a Type la supernova (jMa Sz Asp den! |2010| ; iMa et al.l |2010[ ) . This simulation 
uses a realistic stellar equation of state and a turbulent flame model, and is typical of more 
realistic CASTRO applications. The domain is 5.12 x 10 8 cm on a side, and is covered with 
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512 64 3 grids. There are two levels of factor two refinement, with approximately 1.8% of the 
domain covered by level 2 grids with an effective resolution of 2.5 x 10 5 cm. Figure [T6l is a 
close-up of the center of the domain so that the level 2 grids are more visible. 



9. Summary 

We have described a new Eulerian adaptive mesh code, CASTRO, for solving the mul- 
ticomponent compressible hydrodynamic equations with a general equation of state for as- 
trophysical flows. CASTRO differs from existing codes of its type in that it uses unsplit 
PPM for its hydrodynamic evolution, subcycling in time, and a nested hierarchy of logically- 
rectangular grids. Additional physics includes self-gravitation, nuclear reactions, and radia- 
tion. Radiation will be described in detail in the next paper, Part II, of this series. 

CASTRO is currently being used in simulations of Type la supernovae an d core-collapse 



supernovae; examples of sim ulations done using CASTRO can be found in IJoggerst et al. 



( 120091 ): IWooslev et all (120091) . Further details on the C ASTRO algorithm can be found in 



the CASTRO User Guide flCASTRO User Guide! l2009h . 
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Fig. 1. — Cartoon of AMR grids with two levels of factor 2 refinement. The black grid covers 
the domain with 16 2 cells. Bold lines represent grid boundaries, the different colors represent 
different levels of refinement. The two blue grids are at level 1 and the cells are a factor of 
two finer than those at level 0. The two red grids are at level 2 and the cells are a factor of 
two finer than the level 1 cells. Note that the level 2 grids are properly nested within the 
union of level 1 grids, but there is no direct parent-child connection. 
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Fig. 2. — Weak scaling behavior of the CASTRO code on the jaguarpf machine at the OLCF. 
For the two-level simulation, the number of cells that are advanced in a time step increases 
by a factor of three because of subcycling. To quantify the overall performance, we note that 
for the 64,000 processor case without gravity, the time for a single core to advance one cell 
for one time step is 24.8 /xs. 
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Fig. 3. — Adaptive CASTRO solution vs. analytic solution for Sod's problem run in ID at 
an effective resolution of 128 cells. 
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Fig. 4. — Adaptive CASTRO solutions vs. analytic solution for the double rarefaction prob- 
lem run in ID at effective resolutions of 128, 512 and 2048 cells. 
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Fig. 5. — Adaptive CASTRO solution vs. analytic solution for the strong shock problem run 
in ID at an effective resolution of 128 cells. 
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Fig. 6. — CASTRO solution at t — 0.01s for the spherical Sedov blast wave problem run in 
ID spherical, 2D axisymmetric, and 3D Cartesian coordinates. This was run with a base 
grid with Ax = 0.03125 cm and 3 levels of factor 2 refinement for an effective resolution of 
Ax = .00390625 cm. 
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Fig. 7. — CASTRO solution at t — 0.01s for the spherical Sedov blast wave problem run in 
ID spherical, 2D axisymmetric, and 3D Cartesian coordinates. This was run with a base 
grid with Ax = 0.03125 cm and 4 levels of factor 2 refinement for an effective resolution of 
Ax = .001953125 cm. 
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Fig. 8. — CASTRO solution at t = 0.1s for the cylindrical Sedov blast wave problem run in 
2D Cartesian coordinates. This was run with a base grid with Ax = 0.03125 cm and 3 levels 
of factor 2 refinement for an effective resolution of Ax = .00390625 cm. 
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Fig. 9. — Density in a single-mode Rayleigh- Taylor simulation for a variety of advection 
schemes. Dimensionally-split method results are shown on the top row; unsplit method re- 
sults are shown on the bottom row. We see that the unsplit methods do better at suppressing 
the growth of high-wavenumber instabilities resulting from grid effects. 
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Fig. 10. — Maximum magnitude of velocity vs. time for the stationary star gravity problem. 
At each of the three resolutions, we show the results using a monopole gravity approximation 
and Poisson solve for gravity. We note that in each simulation, the maximum velocity is 
not strictly increasing, leading us to believe that over longer periods of time the velocities 
will remain small. We note that sound speed at the center of the star is approximately 
9.4 x 10 8 cm/s, so at the highest resolution, the peak velocity is less than 1% of the sound 
speed. The solutions in the monopole and Poisson cases match up very well; the discrepancy 
we see at the finest resolution is due to large velocities at the edge of the star, which is 
typically outside the region of interest. 
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Fig. 11. — Magnitude of velocity vs. radius at t — 1 s for the stationary star gravity problem. 
At each of the three resolutions, we show the results using a monopole gravity approximation 
and Poisson solve for gravity. Here, we see clear second order convergence in the max norm, 
and the monopole and Poisson simulations agree best at the highest resolution. 
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Fig. 12. — Radius vs. time for the homologous dust collapse problem in ID, 2D, and 3D 
simulations as compared to the exact solution. In all three cases we see excellent agreement 
with the exact solution. 
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Fig. 13. — Comparison of FLASH (red) and CASTRO (blue) temperature contours for the 
reacting bubble test. Temperature contours at 10 8 , 1.5 x 10 8 , 2 x 10 8 , 2.5 x 10 8 , 3. x 10 8 , 
3.5 x 10 8 , 4. x 10 8 , 4.5 x 10 8 , 5. x 10 8 , 5.5 x 10 8 , 6. x 10 8 , 6.5 x 10 8 , 7. x 10 8 , 7.5 x 10 8 , 
8. x 10 8 K are shown, drawn with alternating solid and dashed lines. The inset shows the 
detail of the middle bubble. We see good agreement between FLASH and CASTRO. 
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Fig. 14. — Comparison of FLASH (red) and CASTRO (blue) 24 Mg mass fraction contours 
for the reacting bubble test. Contours are drawn at values of X of 5 x 10~ 9 , 5 x 10~ 8 , 
5 x 10~ 7 ,5 x 10~ 6 , with alternating solid and dashed lines. The inset shows the detail of 
the middle bubble. As with the temperature, we see good agreement between FLASH and 
CASTRO. 
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Fig. 15. — Here we see a 2D slice of the temperature field from a 3D calculation of a Type 
la supernova with two levels of refinement. There are 512 grids, each containing 64 3 cells, at 
the coarsest level, over 1000 grids at level 1 and over 2000 grids at level 2. Approximately 
1.8% of the domain is at the finest resolution. 
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Fig. 16.- 
grids. 



Here we see a close-up of the previous figure, showing more detail of the level 2 



